Distribution Of Wealth In A Network Model Of The Economy 
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I. INTRODUCTION 
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We show, analytically and numerically, that wealth distribution in the Bouchaud-Mezard network model of 
the economy is described by a three-parameter generalized inverse gamma distribution. In the mean-field limit 
of a network with any two agents linked, it reduces to the inverse gamma distribution. 
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Wealth distribution has become a subject of keen interest in econophysics research. (l| Here, we study a network model of the 
economy proposed by Bouchaud and Mezard (BM).|2| In the mean field (MF) limit of a completely connected network, where 
. any two agents in the network are linked, the model yields the inverse gamma (IGa) stationary wealth distribution. An important 
t— I " feature of the IGa distribution is the power-law (PL) tail. iQ] In the opposite limit of a completely disconnected network, the time- 
dependent part of wealth distribution is lognormal (LN). Both LN and IGa have long history in models of wealth distribution. 
LN is generated by Gibrat's law.|@] IGa, as well as a specific form of GIGa (generalized inverse gamma distribution), were used 
to analyze wealth distribution in ancient Egypt. 

Souma et al. did a numerical study of the BM model and proposed that there may be quite an abrupt transition between IGa 
to LN as a function of the number of connections and the type of connections - regular network or a small-world network. [6] 
C^h . In this paper, we revisit Souma's simulations and compute the p-values of the fitting distributions using Kolmogorov-Smirnov 
test J3|. We argue that the time-dependent LN distribution is a transient - albeit possibly slow, depending on the parameters 
- and concentrate on the stationary solution. We find that for the BM model the latter is a three-parameter GIGa distribution. 
Theoretically, we develop an effective field theory for the BM model of partially connected networks, including regular network 
and a random small-world network IQ] and obtain the Fokker-Planck equation (FPE) for the probability density function (PDF). 
Its stationary solution is a GIGa distribution, with IGa distribution as its limit in the MF regime. 

This paper is organized as follows. In Sec. [TTJ we discuss the effective field theory of the BM model, the corresponding 
\Q • stationary FPE and its GIGa solution. In Sec. [TTTJ we present the results of our numerical simulations. In Sec. IIV1 we summarize 
fN| ' our findings and outline future directions of our work. 
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^ ; II. THEORY 

A. GIGa from Bouchaud-Mezard model 

The BM model reads Q]: 

dWj = V2oWidBi+ £ JijWjdt- £ JjiWjdt, (1) 

jm ' im 

where = means that the stochastic differential equation (SDE) is interpreted in the Stratonovich sense 10,0] and ; = 1,2, ...,N 
with N ^ \ the total number of agents, W, is the wealth of an agent, dBj is an independent Wiener process and a and 7,j are 
constants. Since the BM model may have a wider range of applications - including possibly neural networks - than originally 
intended thus we will study it without applying specific interpretations to W and the model parameters. 
The BM model in (Q} can be rewritten into an Ito SDE H 

dWj = VloWidBi + a 2 Widt + £ JijWjdt - £ JjiWidt. (2) 
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in agreement with Souma et al. @]. Rescaling per W,(f ) = wi{t)e a ', we obtain 



dwi = y/lOWidBi ; + ^ JijWjdt— ^ JpWidt. 



(3) 



It is easily seen that in the large N limit, dwi = V2<yY^=i WidBi ss 0, which implies that the total "wealth" fluctuates around 
a constant value. 

Ultimately, the goal is to determine the PDF P(w,t). Towards this end we notice that there is a discontinuous transition from 
the interacting case /y 7^ to the non-interacting case /y = 0, that is, as soon as the interaction between the agents is turned on, 
the nature of the distribution function is qualitatively changed. Indeed, as follows from eq. (7.8) in J3], P(w,t) does not have a 
stationary limit for/y = and decreases to zero for any finite w when t — s- +°° (while preserving the total "wealth"): 



P(w,t) 



1 



iJntaw 



exp 



1 /logw + a 2 f 

2 V V2t<7 



(4) 



Conversely, in the /y 7^ 0, a stationary solution P(w) = P(w,°°) exists and in what follows we concentrate on its analytical 
derivation while leaving dynamics to numerical investigation. 

The MF limit of a completely connected network was studied in |2j]. Substituting /y = J/N in (01 and extending summation 
on j to each member of the network, we obtain 



dwi = V2awidBi + J(w — Wi)dt , 



where w = N 1 Y,f=i w i is the average of w,. The corresponding FPE is given by 



dP(w) _ d[J{w-w) + o 2 w]P 
dt dw 



dwP 



w- 



(5) 



(6) 



Rescaling via w — > w/w, so that w — 1, we find the normalized stationary IGa solution J2j] 



P(w) 



-2-- 



J+a 2 



(7) 



with a PL tail oc w-P+J/o 2 ), 

For a partially connected network, where each agent is connected with \<n = zN<{N — 1) other agents (0 < z < 1), we 
substitute /y = J/n in (01 and notice that 



£ J ij (Wj-W i )=J(wW-W i ), 
interacting agents: 



(8) 



where w W = n 1 Linteiacting agents: y'(^i) w j I s me average over interacting agents. Observing that h^"' ~ W; when n ~ 1 and 
vjjW ~ 1 when n ~ N, we introduce an effective field theory ansatz: 

->0(y)w 1 -'', (9) 

where 7=1 corresponds to the MF limit and 7 — >• to the minimally connected network. The corresponding FPE becomes 

dwP~ 



dP(w) _ d[J{w-ew l -^) + a 2 w]P + ^ 2 d_ 
dt dw dw 



w- 



(10) 



which has the normalized GIGa solution 



P(w) 





y 2 


r 







try w 



-2-- 



(11) 



with the same PL tail « w ( 2+J l al ) as we find in the MF limit. 6(7) is determined from the normalization condition 

wP(w)dw = 1, 



W : 



(12) 
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and is given by 



6(7) = 



7<t z 



Gamma 



Gamma 




(13) 



It is a monotonic functions between the endpoints 0(1) = 1 and 0(0) = 1 + <J 2 /J respectively. 

We emphasize that (fTTT i describes a GIGa distribution for any y G (0, 1] with the limit of IGa for 7=1. 7 has a finite 
lower cut-off even if only a few connections for each agent are present, while 7=0 corresponds to a completely disconnected 
network which does not have a stationary solution. The latter was explained before but also follows from ( fTQb since the term in 
parentheses (multiplying J) in the r.h.s. is zero in the 7=0 limit. 



III. NUMERICAL SIMULATION OF BOUCHAUD-MEZARD MODEL 



We employ the numerical algorithm described in the Appendix. The time evolution of the distribution function (and its 
parameters) is observed on approach to its stationary limit. Two cases of the network model are considered: Jg] 

1 . In a regular network (RegN), each agent connects with 1 <n=zN < (N — 1 ) nearest neighbors on a circle. In the simplified 
model, we set = J/n in the numerical simulation. 

2. In a random small-world network (RanN), any two agents on a circle have a probability psw to be connected. In the 
simplified model, we set Jjj = J/n in the numerical simulation, where n = pswN. lfl7ll 

Not surprisingly, the main differences between the two networks are as follows: 

• RanN has considerably shorter equilibration time than RegN on approach to a stationary distribution; 

• RanN is better described by EFT than RegN. 
The GIGa PDF is given by 

GIGa(a,j3, 7 ;w) = ^— e W (H) . (14) 
P r («) \wj 

While in theoretical description above the mean is set to unity ( fT2b , which would stipulate j8 = F(a)/T(a — 1/y), numerically 
we fit with a three-parameter GIGa, which allows deviation from w = 1. In our simulation we use a = Vom and J = 0.1. 
Comparison of ( fT~4T > with (fTTT i yields 



ay = l+//a 2 = 3, (15) 

whence in the MF theory limit, y = 1, we have a = 3 and j3 = 2. 

The results of our numerical simulations are presented as follows. Distribution fitting for RegN and RanN are shown respec- 
tively in Figs. Q]and[2] Time evolution of the parameters for RegN and RanN are shown respectively in Figs. [3]and|4] Fig. [5] 
shows the time evolution of the parameters in the EFT and is color coded to be congruent with Fig. @] 

It should be mentioned that we used both the p-value and the log-likelihood to measure which fitting distribution is better, but 
since they yield the same results, in Figs. [3]|5]we present only the former. Also, since the LN fit fails outside initial short time 
scales for parameters used here, we omit it from these plots as well. Figs. [3] and [4] clearly demonstrate that GIGa provides a 
better fit than IGa except on approach to MF limit. It is also clear that the equilibration time for establishment of parameters y 
and ay is much faster for RanN than for RegN. lfl8ll Comparison between Figs. [4]and[5]shows that the EFT describes RanN very 
well. 



IV. DISCUSSION 



We demonstrated that the stationary solution of the partially connected Bouchaud-Mezard network model of the economy 
predicts the generalized inverse gamma distribution of "wealth." In the mean-field limit of a fully connected network, we 
recover the inverse gamma distribution obtained in J2]. lfl9ll Interestingly, the generalized inverse gamma also describes well the 
distribution of human response times lfl5ll . We speculate that for some cases where Pareto distribution is observed it could be, 
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Wealth Wealth 

FIG. 1: Histogram of "wealth" distribution for regular network with z = 0.01 at t = 1 (left) and t = 2500 (right). Red: LN; 
Green: IGa; Blue: GIGa. The fitting of LN is good only for very short times. At t — 1, the p-value of LN is 0.18; of IGa is 0; of 
GIGa is 0. At t = 2500, the p-value of LN is 0; of IGa is 0; of GIGa is 0.05. For t > 2, the p-value of LN is always 0. 
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FIG. 2: Histogram of "wealth" distribution for random small world network with p$w = 0.003 at f = 1 (left) and t = 500 
(right). Red: LN; Green: IGa; Blue: GIGa. The fitting of LN is good only for very short times. At t = 1, the p-value of LN is 
0.51; of IGa is 0; of GIGa is 0. At t = 500, the p-value of LN is 0; of IGa is 0.006; of GIGa is 0.97. For t > 3, the p-value of LN 

is always 0. 



in fact, a power-law tail of a (generalized) inverse gamma distribution. For instance, the distribution function of landslide area 
ind. 

We described partially connected networks using an effective field theory and showed that it describes a randomly connected 
small world network particularly well, including the transitory behavior. Applicability of random versus regular network may 
depend on the circumstance, given the wide relevance of BM model, ranging from psychology to economics. 

An interesting extension of this work would be to study the Watts and Strogatz network model IU2I1 . which has both clustering 
and small-worldness properties. (This model has already been discussed by Souma et al. 10].) The Watts and Strogatz model is 
more realistic than RanN as any real economic network should reflect both clustering and small-worldness nature. 



Appendix A: Numerical simulation method 

A comprehensive list of numerical simulation methods of SDEs is given in 111 411 . In this paper, RanN and EFT are simulated 
by Milstien's method, which has one-order accuracy, and RegN is simulated by the order 1 .5 strong Taylor scheme. Below we 
only present Milstien's method, which can be found in jjl and fjj) . In Milstien's method, the SDE 



is written into the difference equation 



For the BM model, g = \/2awi, hence 



dx = f(x,t)dt + g(x,t)dB, (Al) 



Ax = fAt + gAB + |^[(AW) 2 -Af]. ( A2 ) 
2 ox 



(A3) 

2 owi 



5 



1.0 ; ij^^t^f,^ 



0.8 

« 0.6 
O 

O 0.4 
0.2 
0.0 



Htiv 




6.0 
5.5 
?s 5.0 

2 4.0 



1000 



1.0 
0.8 
0.6 
0.4 



2000 3000 4000 
time 



5000 



2 0.2 
0.0 




1000 2000 3000 4000 5000 
time 




1000 2000 3000 4000 5000 
time 



l.Or 
§ 0.8 

1 0.6 
a 0.4 

2 0.2 
0.0 



1000 2000 3000 4000 5000 
time 



3.5 
3.0 
2.5 




1.0 

I 0.8 

T 0.6 

m 0.4 
O 

O 0.2 
0.0 



1000 2000 3000 4000 5000 
time 

••r»».! 




1000 2000 3000 4000 5000 
time 




1000 2000 3000 4000 5000 
time 



i 

a. 

si 

O 



1.0 
0.8 
0.6 
0.4 



3 °- 2 

0.0 



< C 7 



1000 2000 3000 4000 5000 



time 



FIG. 3: Regular network. First row: y of GIGa (colored), ay of GIGa (colored). Red, green, and blue for z = 0. 1 , 0.01 , 0.003 
respectively. Second, third, and fourth rows (uncolored) are p-values forz = 0.1,0.01,0.003 respectively. 



and 



Ax = fAt + gAB + a 2 Wi[{AB) 2 - At] 



(A4) 



Acknowledgments 

We wish to thank Prof. Dr. Igor Sokolov for many helpful discussions at the early stage of this work. 



[1] V.M. Yakovenko and J.B. Rosser Jr. Colloquium: Statistical mechanics of money, wealth, and income. RMP, 8 1(4): 1703, 2009. 
[2] J.P. Bouchaud and M. Mezard. Wealth condensation in a simple model of economy. Physica A, 282(3):536-545, 2000. 
[3] Y. Fujiwara, W. Souma, H. Aoyama, T. Kaizoji, and M. Aoki. Growth and fluctuations of personal income. Physica A, 321(3):598-604, 
2003. 



6 



[4] H. Aoyama, Y. Fujiwara, and Y. Ikeda. Econophysics And Companies: Statistical Life And Death In Complex Business Networks. 

Cambridge Univ Pr, 2010. 
[5] A. Y. Abul-Magd. Wealth distribution in an ancient egyptian society. PRE, 66(5):057104, 2002. 
[6] W. Souma, Y. Fujiwara, and H. Aoyama. Small-world effects in wealth distribution. ArXiv cond-mat/0108482, 2001. 
[7] Aaron Clauset, Cosma Rohilla Shalizi, M. E. J. Newman. Power-law distributions in empirical data. arXiv.0706.1062, 2007. 
[8] K.Jacobs. Stochastic Processes For Physicists: Understanding Noisy Systems. Cambridge Univ Pr, 2010. 

[9] Weibull.com. The Generalized Gamma Distribution and Reliability Analysis http://www.weibull.com/hotwire/issuel5/hottopicsl5.htm . 
[10] J.F. Lawless and JF Lawless. Statistical models and methods for lifetime data. 1982. 
[1 1] M. Evans, N. Hastings, and B. Peacock. Statistical distributions. 2000. 

[12] D.J. Watts and S.H. Strogatz. Collective dynamics of 'small-world' networks. Nature, 393(6684):440^142, 1998. 
[13] Gavin E. Crooks. The Amoroso Distribution arXiv:1005.3724vl, 2010 

[14] Peter E. Kloeden and Eckhard Platen. Numerical Solution of Stochastic Differential Equations. Springer- Verlag, 1992. 
[15] Tao Ma, R. A. Serota and John G. Holden (unpublished) 

[16] Donald L. Turcotte, Bruce D. Malamud Landslides, forest fires, and earthquakes: examples of self-organized critical behavior Physica 

A, 340:580-589, 2004. 
[17] We use the symbol psw to be consistent with @]. 

[18] The initial condition in our simulations was chosen either w/(f = 0) = 1 or from a narrow Gaussian distribution of w; centered around 1. 

Since in the MFT w<") = w, such choice of initial conditions favors IGa at short times. In Fig[3] one observes slow equilibration of the 

RegN from y = 1 of IGa to y < 1 of the respective GIGa distributions. 
[19] It should be mentioned that the GIGa and IGa are members of a family of distributions related by parametric transformations. (9 UI3I1 

For instance, when a tends to +°° quadratically and y tends to + linearly, GIGa tends to a LN distribution. This limiting circumstance 

corresponds to J /a 2 — > +°°. Thus, by tuning parameters J and a, the BM model may be able to generate a stationary LN distribution. 



7 




400 600 800 1000 
time 




200 400 600 800 1000 
time 



l.Or 

§ 0.8 

1 0.6 
* 0.4 

03 

2 0.2 



0.0 



200 400 600 800 1000 
time 



l.Or 

§ 0.8 

1 0.6 
a 0.4 

2 0.2 
0.0 



l.Or 

§ 0.8 

1 0.6 
a 0.4 

CO 

2 0.2 
0.0 



200 400 600 800 1000 
time 



200 400 600 800 1000 
time 



400 600 800 1000 
time 



l.Orr^ 



1 8 ^^^v^:<^:H 

a 0.4 gp^-tf # v . "jk.¥ 



o.o 



200 400 600 800 1000 
time 

l.O r • . ■ - -a 

8 j-,.:'. " 5-.: « 
0.6 •</*>«.' V * f ' : "\z 



« 0.4 
O 

O 0.2; 
0.0 



200 



400 600 
time 



800 1000 



l.Or 

0.8 
0.6 
0.4 



' 200 400 600 800 1000 
time 



i 

c 
id 
O 

O 



l.Or 

0.8 
0.6 
0.4 
0.2 
0.0 



200 400 600 800 1000 
time 



FIG. 4: Random small world network. First row: 7 of GIGa (colored), ay of GIGa (colored). Red, orange, green, and blue for 
Psw = 0.1,0.003,0.002,0.001 respectively. Second, third, fourth and fifth rows (uncolored) are p-values for 

psw = 0.1,0.003,0.002,0.001 
respectively. 
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FIG. 5: Effective field theory. First row: y of GIGa (colored), ay of GIGa (colored). Orange, green, blue, and purple for 
Yeft = 0.8, 0.6, 0.5, 0.4 respectively. Second, third, fourth and fifth rows (uncolored) are p-values for Jeft = 0.8,0.6,0.5,0.4 

respectively. 



